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Abstract 

In  this  paper  we  describe  a  novel  approach  for  the  solution  of  inviscid  flow  problems  for  subsonic 
compressible  flows.  The  approach  is  based  on  canonical  forms  of  the  equations,  in  which  subsystems 
governed  by  hyperbolic  operators  are  separated  from  those  governed  by  elliptic  ones.  The  discretizations 
used  as  well  as  the  iterative  techniques  for  the  different  subsystems,  are  inherently  different.  Hyperbolic 
parts,  which  describe,  in  general,  propagation  phenomena,  are  discretized  using  upwind  schemes  and 
are  solved  by  marching  techniques.  Elliptic  parts,  which  are  directionally  unbiased,  are  discretized 
using  h-elliptic  central  discretizations,  and  are  solved  by  pointwise  relaxations  together  with  coarse  grid 
acceleration.  The  resulting  discretization  schemes  introduce  artifiw  al  viscosity  only  for  the  hyperbolic 
parts  of  the  system;  thus  a  smaller  total  artificial  viscosity  is  used,  while  the  multigrid  solvers  used  are 
much  more  efficient.  Solutions  of  the  subsonic  compressible  Euler  equations  are  achieved  at  the  same 
efficiency  as  the  full  potential  equation. 
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1  Introduction 


In  the  past  decade  a  substantial  effort  has  been  invested  in  solving  the  Euler  equations,  with  multigrid 
methods  playing  a  central  role.  The  two  major  directions  of  research  in  multigrid  solution  of  the  Euler 
equations  are  the  use  of  coarse  grids  to  accelerate  the  convergence  of  the  fine  grid  relaxations,  and  the  use 
of  defect  correction  eis  an  outer  iteration,  coupled  with  use  of  multigrid  to  solve  for  the  low  order  operator 
involved  [3].  Extensive  research  has  been  conducted  in  both  directions.  Unfortunately,  both  approaches 
can  be  shown  to  have  limited  potential.  Methods  bcised  on  defect  correction  have  /i-dependent  convergence 
rates  for  hyperbolic  equations,  where  h  is  the  mesh  spacing;  other  multigrid  methods  have  p-dependent 
convergence  rates,  where  p  is  the  order  of  the  scheme  involved.  This  unacceptable  situation  motivated  the 
research  outlined  here. 

The  poor  behavior  of  coarse  grid  acceleration  for  hyperbolic  equations,  which  becomes  even  worse  with 
high  order  discretizations,  leads  us  to  conclude  that  coarse  grids  should  not  be  used  to  accelerate  convergence 
for  hyperbolic  problems.  Rather,  the  relaxation  process  alone  should  converge  all  components  of  such 
problems.  This  is  possible,  since  hyperbolic  problems  describe  propagation  phenomena  for  which  marching 
in  the  appropriate  direction  is  a  highly  effective  solver.  For  elliptic  problems,  on  the  other  hand,  local 
relaxation  with  good  smoothing  properties  can  be  greatly  accelerated  by  coarse  grid  correction.  Moreover, 
elliptic  problems  cannot  be  solved  efficiently  by  local  relaxation  alone,  so  that  coarse  grid  acceleration  is 
essential.  In  short,  hyperbolic  equations  do  not  need  coarse  grid  acceleration,  while  elliptic  equations  require 
such  acceleration. 

These  observations  have  motivated  a  study  concerning  the  separation  of  the  hyperbolic  and  elliptic  parts 
in  steady  state  inviscid  flow  calculations.  The  result  is  a  canonical  form  for  the  inviscid  equations,  where  the 
hyperbolic  and  elliptic  parts  reside  in  different  blocks  of  an  upper  triangular  matrix  for  the  system  [5].  These 
canonical  forms  are  the  analog  of  the  decomposition  of  the  time  dependent  one-dimensional  Euler  equations 
into  characteristic  directions  and  Riemann  invariants.  The  insight  gained  by  the  use  of  the  canonical  variables 
enables  one  to  construct  genuinely  multidimensional  schemes  for  these  equations.  It  unifies  the  treatment  of 
the  compressible  subsonic  case  with  the  incompressible  case,  although  these  two  cases  have  been  studied  by 
different  methods  until  now.  Canonical  boundary  conditions  [5]  which  enable  the  proper  numerical  treatment 
of  general  boundary  conditions  are  also  obtained.  Moreover,  these  canonical  forms  enable  the  construction 
of  new  solvers  having  much  better  efficiency  than  existing  solvers. 

The  new  discretization  schemes,  which  are  based  on  these  canonical  forms,  use  upwind  discretization 
only  for  the  hyperbolic  variables,  and  use  a  central  h-ellipiic  discretization  for  the  elliptic  variables.  This 
gives  a  better  representation  for  the  physical  phenomena,  since  elliptic  problems  do  not  have  a  bias  in 
any  spatial  direction,  a  property  that  should  hold  true  for  the  discretization  as  well,  if  possible.  Upwind 
discretization  for  hyperbolic  problems,  on  the  other  hand,  is  compatible  with  the  bias  of  information  flow  in 
the  physical  problem.  The  canonical  form,  therefore,  allows  for  a  better  treatment  of  the  inviscid  equations. 
The  resulting  schemes  are  also  compatible  with  the  uniqueness  properties  of  the  inviscid  equations  under 
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different  geometries  and  boundary  conditions.  In  particular,  the  nonuniqueness  of  solutions  for  exterior  flows 
around  smooth  bodies  is  clearly  evident  with  these  schemes;  only  the  addition  of  a  global  condition,  such 
as  circulation  at  infinity,  insures  uniqueness.  In  existing  schemes  this  issue  is  much  more  obscure  and  there 
seems  to  be  no  direct  analog  of  the  physical  behavior. 

We  discretize  the  elliptic  part  of  the  system  with  staggered  grid  schemes,  which  are  h-elltpUr.  The 
importance  of  h-ellipUc  discretizations  is  that  they  yield  solutions  free  of  spurious  or  weakly  spurious  modes, 
in  the  limit  as  the  Mach  number  goes  to  zero.  The  construction  of  fast  solvers  for  &uch  discretizations  is 
well  understood,  even  in  the  case  of  systems,  and  is  quite  straightforward.  By  contreist,  the  non-staggered 
schemes  commonly  used  would  require  special  techniques  for  the  construction  of  optimally  efficient  multigrid 
solvers  (having  efliciency  comparable  to  that  of  the  full  potential  equation).  Such  solvers  require  multiple 
coarse  grids  on  each  level,  even  for  regular  grids  having  cell  aspect  ratios  near  one  [2]. 

The  elliptic  and  hyperbolic  parts  of  the  equations  are  treated  differently  by  the  solver.  Unlike  existing 
solvers,  which  use  coarse  grids  to  accelerate  the  hyperbolic  part  of  the  system  as  well,  the  new  method 
computes  the  hyperbolic  part  via  relaxations  based  on  the  canonical  forms.  This  involves  marching  in  the 
streamwise  direction  for  the  hyperbolic  quantities,  which  are  the  entropy  s  and  total  enthalpy  H .  The  rest 
of  the  unknowns,  e.g.,  the  velocity  components,  are  relaxed  by  a  Kaczmarz  relaxation  using  preconditioned 
residuals,  yielding  a  smoothing  rate  identical  to  that  for  the  full  potential  equation. 

Numerical  results  are  given  for  a  two  dimensional  flow  around  an  ellipse  and  flow  in  nozzle.  These 
problems  already  include  the  major  difficulties  in  real  problems  and  serve  as  a  good  test  for  the  method 
proposed.  Second  order  schemes  are  used  for  both  cases  and  the  solutions  are  obtained  with  convergence 
rates  similar  to  that  of  the  full  potential  equation  (although  the  work  involved  here  is  larger  accounting  for 
the  multiple  equations). 

2  Canonical  Forms  and  Discretization  Rules 

The  discretization  and  efficient  solution  of  elliptic  systems  of  partial  differential  equations  is  quite  well 
understood.  One  of  the  important  concepts  here  is  h-ellipticity  [1].  It  guarantees  that  the  stability  of  high 
frequencies  for  the  discrete  problem  is  in  correspondence  to  that  of  the  differential  system.  For  the  latter, 
ellipticity  is  defined  in  terms  of  the  symbol  P{,u)  as 

det  P(w)  >  (2.1) 

while  h-ellipticity  is  defined  as 

det  P^(e)  >  C\9\^”',  Ifi*!  <  IT  (2.2) 

Discretizations  which  are  k-elliptic  admit  local  relaxation  methods  with  good  smoothing  properties.  This, 
together  with  efficient  coarse  grid  acceleration  for  smooth  components,  makes  standard  multigrid  methods 
very  efficient  for  such  discretizations.  Although  other  types  of  discretization  also  admit  fast  multigrid  solvers, 
we  restrict  our  focus  to  h-elliptic  discretization  for  elliptic  problems. 
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The  discretization  of  hyperbolic  equations,  which  usually  describe  propagation  phenomena,  can  be  done 
naturally  using  upwind  bicised  schemes.  However  the  application  of  the  above  ideas  to  the  steady  state 
inviscid  incompressible  and  subsonic  compressible  equations  is  not  straightforward,  since  these  equations  are 
neither  elliptic  nor  hyperbolic,  but  rather  mixed  hyperbolic-elliptic.  The  optimal  treatment  of  the  problem 
should  therefore  include  an  identification  of  these  two  parts,  which  have  inherently  different  behavior  and 
call  for  different  numerical  treatment  both  on  the  level  of  the  discretization  and  the  solver.  The  device  for 
this  is  a  canonical  form  of  the  equations,  described  in  detail  in  [5]. 

For  a  two  dimensional  flow,  the  canonical  form  of  the  compressible  Euler  equations  is 


where: 

D\  =  pIc^{{c^  -  u^)Dx  -  uvDy) 

D2  =  -  v^)Dy  -  uvDx)  (2.4) 

Dq  =  vDx  ~  uDy 

In  view  of  these  canonical  forms  we  can  use  the  following  discretization  rule  for  the  invisicd  equations; 

1.  Use  central  (unbiased)  h-elliptic  discretizations  for  elliptic  subsystems 

2.  Use  upwind  biased  schemes  for  hyperbolic  subsystems 

3  Discretization 

Let  a  domain  il  6  IR?  be  divided  into  arbitrary  cells.  Let  the  vertices,  edges,  and  cells  be  V,E,  and  C 
respectively.  The  well  known  Euler  formula 

#l/q.#C-|-#ho/es=  #£-1- 1  (3.1) 

suggests  several  possibilities  for  discretization  of  different  systems  on  structured  and  unstructured  meshes. 
Rewriting  the  Euler  formula  as 

#v  +  # V  -b  #C  +  #C’  +  #/io/es  =  #(7  -b  #7  +  #£  +  1  (3.2) 

one  obtains  the  following  choice  of  discretization.  Let  H  be  associated  with  the  cell  centers,  the  normal 
velocity  components  be  at  the  edges,  and  the  entropy  be  at  the  vertices.  Quantities  other  than  these 
are  calculated  by  well  known  algebraic,  relations  for  the  thermodynamic  quantities  and  by  averaging  for 
the  tangential  velocity  components.  With  each  cell,  we  associate  one  continuity  equation  and  one  energy 
equation,  while  the  two  momentum  equations  are  cissociated  with  each  vertex.  The  following  diagram  is  then 
obtained: 
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The  control  volumes  for  an  unstructured  mesh,  for  the  continuity,  energy,  and  momentum  equations  are 


shown  in  Figure  2. 

The  canonical  form  for  the  compressible  equations  suggests  that  only  the  entropy  and  the  total  enthalpy 
will  be  discretized  using  upwind  biased  schemes,  and  only  in  the  appropriate  terms;  that  is,  only  those  terms  in 
which  a  derivative  in  the  streamwise  direction  occurs  will  be  discretized  with  upwind  bias.  Other  derivatives 
of  these  quantities  will  be  discretized  using  central  differencing.  Let  ei  =  {u/q,v/q),  and  ei  =  {-v/q,u/q). 
Decomposing  the  unit  normal  vectors  to  the  edges  as 


n=  (n  61)61  +  (0  -62)62 


(3.4) 


we  obtain  the  following  discretization  of  the  pressure  terms 


■  ®2)4]‘^si  (3-5) 

16^ 

(3.6) 
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where 

(3,7) 

/  =p(5^//^<,2).  (3.8) 

and  superscript  ”  up”  refers  to  an  upwind  biased  approximations,  while  superscript  ”  c”  refers  to  a  central 
approximation.  The  full  discretization  is  then 

5^PjVj  Ujtfsj  =  0,  (3.9) 

je-r 

^((piVi  •  n,)V(ds(  +  dp  =  0,  (3.10) 

'€7 

X] =0,  (3.11) 

je7 

where  p  is  computed  using  a  symmetric  formula  at  subsonic  points,  and  using  an  upwind  biased  one  for 
supersonic  points.  All  other  quantities  involved  use  central  approximations.  Control  volumes  for  the  different 
equations  are  denoted  by  7  and  7.  Thus  the  only  quantities  that  have  been  upwinded  are  the  entropy  and 
total  enthalpy.  This  results  in  a  reduced  artificial  viscosity. 

4  Multigrid  Algorithm 

The  multigrid  solver,  like  the  discretization,  is  based  on  the  canonical  forms  mentioned  in  section  2.  The 
main  ingredient  differing  from  other  multigrid  methods  is  the  relaxation  scheme.  Other  elements  of  the 
multigrid  method  are  standard  and  will  be  mentioned  only  briefly. 

As  can  be  seen  from  the  canonical  form,  the  hyperbolic  and  elliptic  parts  for  subsonic  flows  are  separated. 
Since  these  subsystems  are  of  very  different  nature,  it  is  unlikely  that  the  same  numerical  process  will  be 
optimal  for  both.  Indeed,  it  can  be  shown  that  coarse  grids  are  inefficient  in  accelerating  certain  smooth 
components  for  hyperbolic  problems.  For  these  components  convergence  rate  is  roughly  (2P  —  l)/2'’  for  a 
p  —  order  method.  The  better  the  scheme,  the  less  coarse  grids  help! 

This  behavior  suggests  that  coarse  grids  are  not  appropriate  for  accelerating  convergence  for  hyperbolic, 
problems.  The  relaxation  should  therefore  converge  all  hyperbolic  quantities  in  the  problem.  While  for 
elliptic  problems  relaxation  cannot  be  efficient  for  smooth  components,  for  hyperbolic  problems  the  situation 
is  different.  Marching  in  the  direction  of  the  physical  flow  is  very  efficient  in  converging  all  hyperbolic 
components  eliminating  the  need  for  coarse  grid  acceleration.  For  elliptic  problems,  on  the  other  hand, 
releixation  techniques  with  good  smoothing  properties,  combined  with  coarse  grid  acceleration,  yield  optimal 
solvers.  The  separation  of  the  different  subsystems  presented  by  the  canonical  form  allows  one  to  construct 
an  optimal  solver  for  the  full  system.  Marching  techniques  will  be  used  for  the  hyperbolic  quantities,  while 
local  relaxation  with  good  smoothing  will  be  used  for  the  elliptic  parts. 
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4.1  Relaxation 


Let  the  residual  of  the  compressible  Euler  equations  be  denoted  by  (R>’,  /f'"',  /f'"  ,  )  and  the  ones  for  the 

canonical  form  by  (r*,r^,r^,r^).  Then  the  following  relation  prevails 


/  1 
—q^A 

\  -W 


0 

— t) 

u 

0 


0 

u 

V 

0 


\ 

/  /?"  \ 

/ 

\  / 

(4.1) 


The  rel^lxation  for  the  Euler  equations  is  employed  as  follows.  The  total  enthalpy  is  relaxed  first,  using 
the  preconditioned  residual  r* 


(4.2) 

Next  the  entropy  is  relaxed  using 

^ +  (»•?)-» (4-3) 

where  and  are  the  diagonal  coefficient  in  the  discrete  version  of  uD^  +  vDy  and  —pT[uDx  +  vDy), 
respectively.  This  is  followed  by  relaxing  the  continuity  and  vorticity  equations.  The  relaxation  for  the 
continuity  equation  is  done  by  keeping  the  values  for  the  density  frozen,  i.e., 

•  Dj  < —  \j  ■  Uj  +  (ri)ypjdsj/  J2{fydskf,  j  e  7,  (4.4) 

and  the  vorticity  equation  is  relaxed  as: 

V(  •  t/  < —  V,  •  t;  +  {rc/{qp)hdsi/  ^{dskf  I  £  7  (4.5) 

*€7 

Note  that  the  preconditioning  of  the  discrete  system  does  not  result  in  an  exact  upper  triangular  form. 
Lower  diagonal  terms  exist,  but  these  are  of  order  O(h^)  and  do  not  affect  the  design  of  the  relaxation  and 
other  numerical  processes. 

The  coarsening  part  of  the  multigrid  method  used  is  standard,  so  its  details  are  omitted.  C^oarse  grids  are 
created  by  combining  neighboring  fine  grid  cells  into  a  coarse  grid  cell.  Linear  interpolation  of  corrections 
and  full  weighting  of  residual  and  functions  are  used  in  an  FMG-FA.S  formulation  [1]. 

5  Numerical  Results 

We  present  here  numerical  results  for  subsonic  cases  of  flow  in  a  nozzle  and  around  an  ellipse.  These  cases 
already  present  most  of  the  difficulties  encountered  in  subsonic  compressible  flows.  For  both  cases  body 
fitted  grid  were  used.  In  the  case  of  the  ellipse  the  grid  extended  to  to  a  distance  of  7  chords. 

Figure  3  and  4  show  results  for  flow  in  a  nozzle.  Mach  lines  as  well  as  relative  errors  in  entropy  are  shown. 
Cases  of  inflow  Mach  number  of  .2  and  .02  are  shown.  Observe  the  small  size  of  the  errors  in  entropy,  as 
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well  as  the  symmetry  of  the  solution.  The  convergence  rale  of  the  multigrid  method  used  is  shown  in  Figure 
5,  in  the  top  two  pictures.  A  reduction  of  eight  orders  of  magnitude  in  10  cycles  is  obtained.  Note  that  the 
rate  of  convergence  is  independent  of  the  Mach  number. 

Figure  5  shows  results  for  a  flow  around  an  ellipse  at  a  10  degrees  angle  of  attack,  with  free  stream 
Mach  number  of  0.1,  with  zero  circulation  at  infinity.  This  is  a  case  of  special  difficulty.  The  exact  values 
for  the  lift  and  drag  coefficient  are  zero.  The  calculated  values  are  Cl  =  2.3  x  10  =  -7.7  X  10-3. 

A  grid  refinement  study  shows  that  the  lift  and  drag  coefficient  approach  zero  with  mesh  refinement.  The 
convergence  history  for  this  problem  is  shown  in  the  bottom  picture  of  figure  6.  The  full  FMCl  history 
(including  coarse  levels)  is  shown.  The  convergence  rate  in  this  Ccise  is  slower  than  for  the  nozzle  flow  cases 
but  still  very  fast  compared  with  existing  methods. 

6  Conclusion 

A  new  approach  for  the  discretization  and  the  solution  of  the  Euler  equations  have  been  presented.  It  is 
based  on  a  canonical  form  of  the  equations  which  allows  for  discretization  which  involve  upwind  biased 
discretization  only  for  the  physically  biased  quantities,  that  is,  entropy  and  total  enthalpy.  The  multigrid 
method  used  for  that  discretization  shows  essentially  optimal  convergence  rates  for  the  second  order  schemes 
used.  Moreover,  nonlifting  solutions  around  smooth  bodies  at  angle  of  attack  are  obtained.  Unlike  most 
other  methods,  the  performance  of  the  method  does  not  degrade  as  the  Mach  number  approaches  zero. 
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Figure  5:  Residulas  History 
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